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ABSTRACT 

We examine the effects of self-gravity and magnetic fields on supersonic turbulence in isothermal 
molecular clouds with high resolution simulations and adaptive mesh refinement. These simulations 
use large root grids (512^) to capture turbulence and four levels of refinement to capture high density, 
for an effective resolution of 8, 196^. Three Mach 9 simulations are performed, two super- Alfvenic 
and one trans- Alfvenic. We find that gravity splits the clouds into two populations, one low density 
turbulent state and one high density collapsing state. The low density state exhibits properties similar 
to non-self- gravitating in this regime, and we examine the effects of varied magnetic field strength on 
statistical properties: the density probability distribution function is approximately lognormal; veloc- 
ity power spectral slopes decrease with field strength; alignment between velocity and magnetic field 
increases with field; the magnetic field probability distribution can be fit to a stretched exponential. 
The high density state is characterized by self-similar spheres; the density PDF is a power-law; collapse 
rate decreases with increasing mean field; density power spectra have positive slopes, P(p, k) oc k] 
thermal-to-magnetic pressure ratios are unity for all simulations; dynamic-to-magnetic pressure ratios 
are larger than unity for all simulations; magnetic field distribution is a power-law. The high Alfven 
Mach numbers in collapsing regions explain recent observations of magnetic influence decreasing with 
density. We also find that the high density state is found in filaments formed by converging flows, con- 
sistent with recent Herschel observations. Possible modifications to existing star formation theories 
are explored. 

Subject headings: methods: numerical — AMR, MHD 



1. INTRODUCTION 

Star formation is one of the most important outstand- 
ing problems in astronomy and astrophysics. Over the 
last six decades, the theory of star formation has gone 
through a major paradigm shift. Early work (Mestel fc| 
|Spitzer||l956[ [Mouschovias 1 1 1 9 76 [ ) focused on the domi- 
nance ot magnetic fields as the primary physical agent 
in star formation, suppressing collapse to support the 
perceived long lifetime of molecul ar clouds. La ter work 
jLarson 1981 ; Elmegreen 1993; Pad oan fc Nordlun d 20021 
[Krumholz McK ee 2005; Padoa rTiriN ordlund] |20lf[ ) 
shifted the locus from magnetically dominated collapse 
to turbulence dominated col lapse. In a swing in the other 
direction, new observations (iLi et al.|2009| |Crutcher et al 



2010[ |Heyer fc Brunt|[2QT2{ [I'argon et al.|[201l[ ) have in- 



dicated that at certain size and density scales, magnetic 
fields dominate, but at smaller scales the field importance 
is reduced. 

This paradigm shift towards turbulence has been made 
possible in large part due to the ever increasing capabil- 
ity of computers and magnetohydrodynamic (MHD) al- 
gorithms, which allow increasingly accurate simulations 
of MHD turbulence (Kritsuk et al.| |2011c^ . Recent al- 
gorithm progress has been made by includir 



ling magnetic 
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fields in high dynamic range codes, mos t notably adap- 



tive mesh r efinement (AMR) (Balsara 200l[ [Fromang 

d 



et al. 1 2 00 6} Collins et al. 2011) and smo othed particle 
hydroaAmam flFrice fc Monaghan 2004 _ Dolag fc Sta- 



Gaburov "STTTitadori 20TT^ 



LTaffents m star i 



The three most important physicalagents m star for- 
mation are gravity, turbulence, and magnetic fields. A 
considerable amount of work has gone into any pair of 
these, but relatively little study has combin ed a l l thre e 
with with high resolution methods. [ Mouschovias (1976), 
Scott fc Black| (|1980D , and IGahi fc Siu, ( ,1993] niaA^elUid- 



ied the combined et ' ects ol magnetic helds a nd gravity. 



Goldreich fc Sridhar] ( |1995D ; |Cho fc Lazarian" (2003 ) and 
Krit suk et al. (2009^ have studied the effects of magnetic 
fields on turbulence, both compressible and incompress- 
ible. Self-gravitating turbulence has b een studied by sev- 



eral au thors, notably Klessenl ( 2000 ) and Kritsuk et al. 
('2011a), finding enhanced high density material relative 
to the pure turbulence results. 

Co mbinin g all three physi cal mechanisnis 



Price fc 



Bate| (|2008[ ), |Federrath et al.| ( [2011b ), and Collins et al 
( |20ll ) have done simulations ol sell- gravitating, magne- 
tized turbulence with high dynamic range methods. The 
primary dif ference in setup between the first of t hose 
works (Price fc Bate 20081 'Federrath e t al.||2011b] ) and 
the work presented here are initial and boundary condi- 
tions: in those works, the initial conditions were isolated 
uniform density spheres with random velocity perturba- 
tions; in the work presented here, we begin with fully 
developed turbulence in a periodic domain. Of course, 
both situations are idealized, with the full nature of star 
formation potentially dependent on the molecular cloud 



2 



formation process as well. 

In this work, we present three high resolution, high 
dynamic range simulations of supersonic, super- Alfvenic 
and trans- Alfvenic turbulence with self gravity. We find 
that gravity breaks the cloud into two distinct states, one 
low density turbulent state and one high density collaps- 
ing state. We discuss several of the dominant statistics 
used in describing supersonic turbulence, the impact that 
the newly formed self-gravitating state has on them, and 
the effect of magnetic fields. In Section [2] we discuss the 
numerical algorithms, initial conditions, and simulation 
parameters. In Section [S] we discuss the density proba- 
bility distribution function (PDF), V{p); followed by the 
density power spectra P{p^k) in Section |4j the distri- 
bution of energy in Section |5] the magnetic probability 
distribution function in Section j6j velocity and magnetic 
power spectra in Section[7| and rmally the relative impor- 
tance of compressible and solenoidal modes in Section [Sj 
We discuss possible implications on star formation, and 
compare with recent observations, in Section |9] We sum- 
marize our findings in Section 10 



2. NUMERICAL METHOD, SIMULATIONS, ANALYSIS 

For the data presented in this paper, we solve the 
ideal MHD equations with self-gravity using the adaptive 



mesh refinement (A MR) code Enzo (Bryan et al. 1995 



O'Shea et aL]|2QQ4[ ) extended to MHD by [Coll ins et al 
( |2U1U|) . This code uses the AMR algorithms d evelop e( 
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by [Berger ^ Cole lla ( 1989p andl Balsara] ( [2QQT| , the hy- 
per ooIiF^oIveF^iLie^ ([2008'), the is otherm al HLLD 
Riemann solve r developea by' Mignonel (2007), and the 
CT method of [ Gardiner fcStone. ( ,200g| . 

We select the Mach number, Al, virial parameter, o^vir, 
and mean thermal-to-magnetic pressure ratio, Po as 



M 



/3o 



^0 



0.2,2,20, 



(1) 



(2) 



(3) 



where I'rms is the r.m.s. velocity fiuctuation, Cg is the 
sound speed, po is the mean density, Lq is the size of the 
box, and Bq is the mean magnetic field. 
These can be scaled to physical clouds as 



tff = l.lnH,V'Myr 

Lo = 4.6cs,2nj^;^3^^pc 

^-1 



= 1.8cs,2km s 



^rms 

M = 5900cs,2nH,V^ ^v/i0 



Bo = (13,4.4,1.3K2nJ/^AiG, 



(4) 

(5) 
(6) 

(7) 

(8) 



where Cs,2 = 0.2km s~^ and nH,3 = uh /{'^000cm~^) are 
the sound speed and hydrogen number density, respec- 
tively, and we have used a mean molecular weight of 2.3 
amu per particle. 
This definition of 

Q^vir IS exact for a uniform density 
sphere in isolation, and we have used it here for consis- 
tency with other works in the literature. In reality, the 
actual importance of gravity relative to kinetic energy 



difficult to ascertain for a turbulent box with periodic 
boundaries due to the infinite nature of the box and the 
intermittent nature of dense structures. The value of 
used here is perhaps somewhat lower th an the aver- 



age value from obs erved molecular clouds (Heyer et al. 



2009 



Dobbs et al.^ 2011), but not outside the observed 
parameter range. Furthermore, the virial parameter is 
potentially scale-dependant, with larger clouds being on 
the average more gravit ationally bound t han smaller ones 
(Heyer et al.||2001( Goodman et al. '2009), so it is possi- 
ble that these results apply better to a subset of large, 
gravitationally bound molecular clouds. 

The initial conditions for this simulation were gener- 
ated by a suite of unigrid s imulations using the PPML 
code fUstyugo v et al.||2009 ) without self-gravity. Cubes 
with 1024^ zones, with initially uniform density and mag- 
netic fields, were driven using a solenoidal driving pat- 
tern. Power in the driving w as between wav enumbers 
k/kmin = 1, 2, and driven as inlMac Low|(|1999|) to main- 
tain our target Mach number. 13" 
eral dynamical times. 



rivmg contmued for sev- 



^dyn — 



in 



2Vr 



1.2njjV^Myr. 



(9) 



until a statistically relaxed state was reache d. Results of 
the tu rbulent boxes were first presented in Kritsuk et al. 
(2009), see Kritsuk et al. (2012, in preparation) for more 
details. 

The simulations were then restarted using Enzo with 
self- gravity and AMR. A root grid of 512^ and 4 levels 
of refinement by a factor of 2 were used, with refinement 
such that the local Jeans length Lj = ^/c^7^/Gp is re- 
solved by at least 16 zones. This gives an effective linear 
resolution of 8, 192. The simulations were run for 0.6tff. 
Figure [l] shows projections through the volume at t = 0.0 
(left column) and t = 0.6tff (right column) for each /3o 
(top to bottom, Po = 0.2, 2, and 20). In the left column, 
one can see the filamentary structures associated with 
supersonic turbulence. The right column shows the dis- 
tribution of high density collapsing cores superimposed 
on the turbulent state. 

It should be noted that the two solvers used for this 
simulation have different dissipation properties. The 
solver used for the initial conditions (PPML) employs a 
third order spatial reconstruction, while the solver used 
in the self-gravitating portion was only second order spa- 
tially. This change in solver was due to the fact that 
AMR as employed in Enzo has not yet been extended 
to include the increased algorithmic complexities of the 
higher quality PPML algorithm. Details about t he differ- 
ences injiumerical dissipation can be found in Kritsuk 
|et al.| ( |2011c| ). In that work, the solver used in these 
AMK runs is referred to as LL-MHD. The effects of this 
solver change can be seen most prominently as a minor 
loss of dynamic range in the velocity power spectra, but 
the statistical properties are otherwise the same between 
the two solvers. 

In order to demonstrate the effects of gravity in the rest 
of the paper, we typically present the simulations at two 
fiducial times, t = O.ltff and 0.6tff. The first snapshot 
was taken at t = O.ltff , which is sufficient to remove the 
effects of the solver transition, but it is early enough to 
not show any effects of gravity. Unless otherwise noted 
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(e.g., velocity power spectra. Figure 15) the statistics at 
t = O.ltff are identical to those at t = O.Otff, which cor- 
responds to the moment at which gravity was turned on. 
Due to the short timescale relevant for the high density 
gas, we average several snapshots around each of the two 
fiducial snapshots, in a range of ±0.05tff , to reduce sta- 
tistical noise. This short-term time averaging is done in 
all plots unless otherwise noted. 

The final time, tgim = 0.6tff, corresponds to O.S^dyn- 
We can estimate the scale at which the turbulence can 
be considered relaxed by using the structure function 
scaling, 6v^ oc ^^p, and computing the scale at which 
the number of turnovers at the end of our simulation is 
greater than some number. A/" = tsim/^£, were t£ = i/Sv£ 
is the turnover time at length £. One finds 



p/(p-Cp) 



l: 



Cp/(P-Cp) 



(10) 



where = 10 is the velocity at the outer scale, Lq = 1 
is the size of the box, and JV is the number of crossings 
at a given scale. For the thi rd order structure fu nction, 
p = 3, in a supersonic fiow, |Kritsuk et al.| ( |2QQ7 ) found 
that = 1-3. For a single turnover time, TV" = 1, we 
find is = 0.09, which corresponds to k/kmin = 11. The 
short time averaging window, ±0.05tff, corresponds to 
k/kmin = 875, which is resolved by all refined regions. 

Power spectra in this simulation were computed only 
using the root grid data, at 512^. Due to the incom- 
plete filling of /c-space, power spectra that also include 
the refined regions would require data interpolation. Due 
to the fact that the volume filling fraction of refined re- 
gions is quite low (see Sectionjsl) , and the fact that power 
spectra are volume weightea quantities, spectra using 
anything but the root grid data would be dominated by 
interpolated data that do not contain much useful infor- 
mation. 

All analysis has been perfor med with the AMR analysis 



2011 ). 



3. DENSITY PDF 

Figure |2] shows the density PDF, V{p)^ for each sim- 
ulation (po = 0.2 in red, /3o = 2 in green, and /3o = 20 
in blue) at t = O.ltff (solid lines) and t = 0.6tff (dotted 
lines). This color scheme will be used throughout the 
paper. The figure additionally shows two grey lines that 
divide the gas into three sections; low density turbulent 
gas (left section), high density self-gravitating gas (cen- 
ter section) , and very high density gas that is numerically 
unresolved (right section) . The first two states are of the 
greatest interest to us here. The last state is interest- 
ing qualitatively, but we cannot make any quantitative 
measurements of the gas here. The transition density be- 
tween turbulent and collapsing states, pt^ is taken where 
the power-law begins to transition from lognormal, at 
Pt = 10. As there is likely shock-compressed gas above 
Pt that is not self-gravitating, pt is not meant to be used 
as a phase boundary between the two states. A more 
complete set of criteria for the transition between turbu- 
lent and collapsing gas is currently under investigation. 
The second division is taken at the highest density that 
is still considered resolved by our refinement criterion on 
the finest level, as discussed in Section |2j This gives 
Pr = 6347. At this density there is also a change in the 



power-law slope. At very high density, corresponding to 
very small scale, inaccuracies in the angular momentum 
transport become dominant and the gas cannot collapse 
fully, leading to excess mass and decreased fragmenta- 
tion. This can be somewhat addressed by incorporating 
sink particles, but as of this writing no satisfactory pre- 
scription of sink particles with magnetic fields has been 
developed. 

In the following, we will identify the turbulent state 
as those features that belong to either low density gas 
(below Pt) oi" gas that whose statistical properties are 
relatively unchanged over the course of our simulation. 
Collapsing gas is identified by high density (between pt 
and pr) and short time variation. 

3.1. Density PDF in the Turbulent State 

One of the most robust properties of isothermal super- 
sonic turbulence is the lognorma l distribut ion of densities 



dBlaisdell et al. 



1993 



zquez- l^jemadenil 1994 
let al.l I1997albt y^calo et a l.| 1^981 Passot & ^ vazquez- 



Padoan 



Vazquez- Semadeni 
[TD98| Passo 

Semadeni 1998; Nordlund & Padoan 1999; Kritsuk et al. 



[2007nFederrath^ al.||2008b| |Price| r2012). Several prop- 
erties of star formation have been predicted using the 
lognormal distribution funct ion, including the initial 
mass function (IMF) of stars ( |Padoan fc Nordlund|20Q2 



Padoan et al. 2007), brown dwarf frequency ( Padoan 



Nordlund 2004) and the star formation ra te (jKrumholz 
& McKee 2005; |Padoan & Nordlund||2011|). 

For compressible turbulence without gravity or mag- 
netic fields, the density PDF, V{p)^ can be shown to be 
a lognormal of the form 



V{p)dlnp: 



1 



V2 



exp 



{Inp- pf 
2a2 



dlnp, (11) 



where /i = —a /2 is the mea n of Inp (Blaisdell et al. 



1993] Vazquez- Semadeni||19 94[) . The variance and Mach 
ber, Al, are related by 



(12) 



The parameter b has been determined numerically t o lie 
between 0.3 and 0.4 (Padoan et al. 1997b; Federrat h et al. 
2008b; Kritsuk et al. 2007; Be etz et a l. 2008; Kritsuk 
,et al.,2010a, ,Federrath et al.,201Qt . The value b has been 
shown t Q"approach unity in simulation s with compressive 
forcing ( [Federrath et al.||2008b [ 



2010). 



The presence of the Lorentz terms in the momentum 
equation breaks the invariance of the equations relative 
to the mean density. Thus the sequence of shocks that 
determine the density of a parcel of gas are no longer 
independent multiplicative events, as they are in hydro- 
dynamic turbulence. For this reason there is no a pri- 
ori expectation of a lognormal density distribution in a 
magnetized system. However, Ostriker et al. (2001) and 
Lemaster & Stone] ( |2008[ ) have shown that V [p) is ap- 
proximately lognornaal, with properties weakly depend- 
ings; on mean field s t rength . For driven MHD turbulence, 
Lemaster & Stone (2008) find that, for densities within 
107o of the peak density, 

(JLsos = VI - 0.72 In (l + 0.5Af2)^ 0.201, (13) 

and weakly decreases with field strength. Kritsuk et al. 
(2012, in preparation) performed high resolution simu- 
lations of statistically stationary magnetized turbulence 



4 




Figure 1. Column density (E) for the two primary snapshots (t = O.Offf, left column, and t = O.Gfff, right column), for all three 
simulations (from top to bottom, (3o = 0.2, 2, 20, respectively). Both figures show the filamentary structure expected from both turbulence 
models and observations. The right column shows the addition of high density collapsing gas, and somewhat enhanced filamentary nature. 
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and showed that the presence of a magnetic field alters 
the low density gas, making a lognormal description less 
appropriate. They did find that the high density wing 
of the PDF is still well approximated by a lognormal for 
the two weak field runs, /3o = 2 and = 20. 

Figure [2] shows the resemblance to a lognormal present 
in our simulations. In Table [T] we show /i and a found 
from fits to the average PDF for three snapshots from 
each of our simulations. Fits were performed for p G 
[5 X 10~^, 10], a range chosen to exclude material above 
pt. We find a weak sensitivity of p and a with both time 
and /3o, with p decreasing with time and increasing with 
Pq. Our short averaging time, relative to a dynamical 
time, means that these can be viewed only as snapshots 
in time, not as robust statistical averages. 

3.2. V{p) Collapsing State 

The PDF of the collapsing material forms a power-law, 
V{p) oc p^, for densities above pt. This was fi r st pre - 
sented in hydrodynamic simulations by [Klessen (2000), 
though the resolution of those simulations was too low to 
measure a significant power- l aw. The first measurement 
of the slope was done by ^Slyz et al. ([2005 ) , who found 
m = — 1.5. This slope was also seen by Federrath et alj 
(|2008al) and iVazquez-Semadeni et al.j (2008) . ,Kritsuk| 
et al.| ( [201 lap used a very high resolution AMR simula- 



tion to measure a slope of —1.67 at intermediate to high 
densities, and —1.5 at high density. They also provided 
and explanation of this power-law by comparing V{p) to 
that of self-similar collapse, wherein p oc r^. The three 
solutions they discussed are the pressure free collapse 
dPenstonl Il96 9[ PF), Larson - Penston supersonic infall 
fLarson 1969i"|Penston 1969[ LP), and expansion wave 



enston 



LP) 



from inside out collapse (" Shu||1977| l EW). Values of n for 
various models, and the implied value of m, are shown 
m Table |2l This table summar izes all semi-analy tic and 
numerical results in this paper. ICollins et al.l (12011) mea- 
sured this exponent for super- Allvenic MHD turbulence, 
and found a value of m=— 1.64 for p G [10,1000]. For 
the simulations in this work, we find m = —1.80, —1.78, 
and -1.65 for (3o = 0.2, 2, and 20. The (3o = 20 case 
presented here is similar in bot h physical para meters and 

and in the 
while 



measured slope to that in Coll ins et al 
gas dynamic simulations of j^Kritsuk et al 



. pllaD 

the stronger field simulations show steeper values. The 
values found here are most consistent with the pressure- 
free value of m = —1.75. 

It is clear from Figure [2] that the slope is a fu nction of 
time in these clouds. This was also discussed by |Kritsuk 



et al. (2011a). Figure [3 shows the time evolution of the 
power- law exponent rnloT the last few snapshots, those 
that contributed to the average shown by the dashed 
curve in Figure |2j The two low field cases, /3o = 2 and 20, 
seem to have converged, with power-law index m in the 
Po = 2 case slightly lower than the Po = 20 simulation. 
The Po = 0.2 simulation has clearly not converged, and 
is still increasing at the end of the simulation. Given the 
agreement of m ^ —1.64 between our 00 = 20 simulation 



and previous simulations (Collins et al. 2011 Kritsuk 



Id p: 

et al.|2011a|) and observations, we feel confident that this 



is a robust result in the hydrodynamic limit. 

As we will discuss in Sections |8] and [9} three dimen- 
sional compressions are suppressed with increasing mag- 
netic field strength, as the flow is forced along magnetic 




Figure 2. Density PDF, V(p), for each of our three-simulations, 
/3o = 0.2, 2, and 20 (red, green, and blue respectively), at two 
snapshots, t = 0.1, 0.6tff (solid and dotted lines, respectively). The 
vertical grey lines separate the low density turbulent state (left 
section), high density collapsing state (center section) and very 
high density unresolved gas (right section). 
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Figure 3. The evolution of the density PDF power-law exponent, 
m, for the /Sq = 0.2, 2, and 20 simulations (red, green, and blue 
respectively). The value of this exponent seems to be increasing 
for the /3o = 0.2 simulation, but possibly reaching a constant value 
for the other two. 



fleld lines. This reduction in compressibility is likely the 
reason for the reduced value of m in /3o = 2, and delayed 
convergence in the /3o = 0.2 run. 

The cumulative mass fraction above some critical den- 
sity Per, deflned as 



M{pcr) = / pV{p)dp, 

J Per 



has been used in a number of theories of star forma- 
tion, as we will discuss in Section [9Jl ln order to quan- 
tify the effects of mean fleld. Figure fi] shows M(pt) for 
each run relative to the /3o = 0.2 simulation. This shows 
an increase in the collapsed mass as a function of fleld 
strength. The behavior of this mass fraction with time. 
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Figure 4. Cumulative mass relative to /3o = 0.2 case for the 
last snapshot, t = 0.6fff. Vertical lines separate turbulent (left) 
collapsing (center) and unresolved (right) gas. 



Table 1 

Lognormal Fit Parameters 





a 




time^~~~~~~~--^^ 


0.2 


2.0 


20 


0.2 


2.0 


20 


0.1 


1.2 


1.3 


1.3 


-0.86 


-0.84 


-0.80 


0.3 


1.2 


1.4 


1.4 


-0.97 


-1.00 


-1.1 


0.6 


1.2 


1.5 


1.5 


-1.1 


-1.2 


-1.4 



where p is the Fourier transform of p, and the star de- 
notes its conjugate. Figure [6] shows the density power 
spectra, P(p, /c), for ah three simulations {/3o = 0.2, 2, 
and 20 in red, green, and blue, respectively) at two snap- 
shots, t = O.ltff and t = 0.6tff (solid, dotted). In Figure 
[7|we show the power spectra for In p, for three snapshots 
of the f3o = 0.2 simulation. Figure [S] shows the column 
density power spectra, P(Yj^k). In Figures [6| and [sj the 
birth of the high density collapsing state carioe seen in a 
dramatic change in the behavior of of P(p, /c), transition- 
ing to a positive slope. This behavior is conspicuously 
absent in the power spectrum of the logarithm of density, 
as we will discuss in the following sections. 




Figure 5. Cumulative mass for p > 10 vs. time for (3q = 0.2,2, 
and 20 (red, green, and blue, resepctively). This shows the rate at 
which gas enters the collapsing state, and is a decreasing function 
of ^0. 



M(p, t), can be used to explain the slower convergence of 
m with the more magnetized simulations. Figure [5] shows 
M{pcr = 10, t), which shows that the rate at which mate- 
rial enters the high density state is a decreasing function 
of Po- This can be used a s a p roxy for star formation, as 
we will discuss in Section [OTH and while the exact rate is 
a function of the critical density pc, the increase of rate 
with Po is not. 

4. DENSITY POWER SPECTRA 

One of the more striking results from this study is the 
behavior of the density power spectrum, P(p, /c), under 
the effects of self-gravity. Here we have defined P(p, k) 
as 



(14) 



4.1. Density Power Spectra in the Turbulent State 

The turbulent initial conditions can be seen in the early 
snapshot (solid line) in Figure [gJ For weakly subsonic 
(but still compressible) isothermal turbulence, one ex- 
pects the density to follow the pressure fluctuations, and 
the density p ower spectrum, P ( p,k) oc fe^, sho uld scale 
as C = -7/3 dBayly et al.||1992| [Kritsuk et al.||2007| ). It 
has been seen that this power spectrum flattens wit h in- 
cre asing Mach number. For trans-sonic turbulence, |Kim| 
|fc Ryu| ( |2005|) measured ( = —1.7 at a Mach number 
\M = 1.2. For sup ersonic turbulence, with = 6, 
Kritsuk et al. (2007) measure ( = —0.9 for a simulation 
with 512^ zones, and a somewhat steeper spectrum of 
( = —1.07 for 1024^. For the early snapshot, we measure 
significantly more shallow spectra, ( = —0.42, —0.58, and 
—0.62 for our /3o = 0.2, 2, and 20 simulations, respec- 
tively. This fit was done for k/kmin ^ [2, 30]. The flatter 
nature of these spectra, relative to the [Kritsuk et al.| 
(2007) work, is likely due to the higher Mach number 
employed in our simulations. For Burger's equation with 
vanishing pressure term, which c ould be analogous to 
increasingly supersonic turbulence, 'Saic hev fc: W oyczyn-| 
(1996) predict that = 0, a trend that is consistent 
with ( " here being shallower than that found in [Kritsuk] 
et al.[ ( [2007[ ). It should be noted that this comparison 



to Burger 's turbulence is not meant to imply that su- 
personic turbulence is similar to turbulence in Burger's 
equation; the lack of vorticity and strong intermittency 
in Burger's equation makes the two systems only super- 
ficially similar. A second potential cause of the change 
in slope is due to the fact that the solver employed here 



is more diffusive that PPM used in Kritsuk et al. (2007 



2011c), which may contribute to the change in slope. I'he 



decrease in slope with increasing magnetic field strength 
is consistent with the decrease in slope in the velocity 
power spectrum, as will be discussed i n Sect ion [7[ 
As discussed in Beresnyak et al. (2005), P[p^k) is 



strongly influenced by the rare, high density peaks. In 
that work, the authors noticed that the power-spectrum 
of the logarithm of density was significantly steeper than 
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P(p, /c), as the logarithm effectively filters out high den- 
sity, rare peaks. This is seen in our simulations as 
well. The behavior of P(ln p, k) shows very little evo- 
lution with time, as the collapsing state is effectively fil- 
tered out, which allows us to recover the turbulent state 
through the entire simulation. This can be seen in Fig- 
ure [7| for the f3o = 0.2 simulation. Slopes get somewhat 
steeper with increasing mean field, with /3o = 0.2 having 
a slope of —1.0, and /3o = 2 and 20 having slopes of —1.4 
and —1.53, respectively (neither shown here). 

In order to make contact with observable quantities, we 
show the power spectra for column density S, P(5], k) ex 
/c^, in Figure [sj For the range k/kmin ^ [2,30] at early 
times, we fina A = -1.39,-1.61, and —1.53 for (Sq = 
0.2, 2, and 20, respectively, roughly consistent wit h the 
addition of —1 to C from the integration. P adoan et al. 
([2004| ) measured P(I],/c), for three star- forming clouds, 
Perseus, Taurus, and the Rosetta nebula. They found 
significantly steeper spectra than we do here, ~ —2.8 for 
all three clouds. They also performed synthetic obser- 
vations through two simulations with = 10, a super- 
Alfvenic simulation and an equipartition model, and find 
A = —2.25 for the equipartition model and —2.7 for the 
super-Alfvenic model. The difference in slope between 
the slopes they found and the slopes in our simulation is 
most likely due to the observations and radiative trans- 
fer models missing the high density material due to the 
limited dynamic range, as ^^CO freezes onto dust grains 



at densities above 
discussed in 



10^ 



Beresnyak et al. 



(IBac mann et al. || 2002|) . 
( |2UU5| and shown in Fi, 



As 

gure 

[7| the slope of the density power spectra are quite sen- 
sitive to the rare high density material, so even a slight 
decrease in the dynamic range of the observations will 
cause a steepening in the spectra. To properly compare, 
we will need to perform similar synthetic observations of 
our simulations. 

4.2. Density Power Spectra in the Collapsing State 

The most prominent effect of gravity is the increase in 
slope of P(p, k) with time, as seen in Figure |6] For self- 
similar spheres with p (x r'^, one finds ( = — 2(n + 1). 
Thus, positive slopes will be seen in P(p, k) for any 
n < —1. Expected values from different self-similar mod- 
els are shown in Table |2] Our measured values of ( for 
^/^min e [10, 200] are 0.86, 1.12, 1.2 for ^^0=0. 2, 2, 20, re- 
spectively, at t = 0.6tff. Direct comparison between the 
measured value of ( and those predicted from the self- 
similar collapse models should be handled carefully for 
two reasons, both stemming from to the volume weighted 
nature of power spectra: first, he measured value con- 
tains contributions from both the collapsing state and 
the turbulent state, which will tend to decrease ( from 
the values expected from pure self-similar spheres; sec- 
ond, the power spectrum contains contributions from un- 
resolved gas at very high densities. 

The effect of increasing mean magnetic field on P(p, k) 
in the collapsing state is to decrease the amount of power 
at all scales, increase the wavenumber at which the slope 
becomes positive, and slightly decreasing the slope. This 
is consistent with the increased compressibility and in- 
creased rate of collapse found in the more weakly mag- 
netized simulations. 

The column density power spectral slope. A, becomes 
nearly flat for the cohapsing case, with A = —0.11,0.15, 




^/^min 



Figure 6. Density power spectra P(p, k) for all three simula- 
tions (/3o=0.2, 2, 20 colored red, green, blue) and two snapshots, 
t = 0.1,0.6tff (solid, dotted). The increasing nature in the later 
snapshot is expected from self-similar spheres. 

and 0.25 for /3o = 0.2, 2, and 20, respectively. This differs 
greatly from what is observed, owing to the fact that the 
material causing the flat spectra is extremely high den- 
sity, with p > {10^ - 10^)po ^ (10^ - 10^°)cm-^, wherein 
typical observational tracers such as CO are frozen onto 
grains, and an extremely small volume filling fraction, re- 
quiring high resolution observations. High density trac- 
ers, such as NH3 or deuterated species suc h as H2D^, 
will be necessary to obser ve this signal (Walmsley et al. 



2004 di Francesco et al. 2007|. Synthetic observations 



of our data, as well as high resolution observations in a 
broad range of chemical tracers, will be required to fur- 
ther reconcile this discrepancy. Again, increased power 
with increased /3o is consistent with increased compress- 
ibility of the gas. 

5. ENERGY RATIOS 

The balance of energies sheds important light on the 
physical processes at play in these clouds. Figure [9] shows 
thermal-to-magnetic pressure ratio, /3th = Sttc^p/B'^^ 
versus density p; and Figure [lO] shows dynamic-to- 
magnetic pressure ratio, ^dyn = Snpv'^/B'^. In both fig- 
ures, the left column is taken at t = O.ltff , the right from 
t = 0.6tff , and top to bottom show Po = 0.2, 2, and 20, 
respectively. Both figures are colored by mass fraction 

5.1. Energy Ratios in the Turbulent State 

In the turbulent state, Pth shows some interesting vari- 
ation with mean field. The scatter in /3th increases with 
increasing /3o, and the mean slopes in the Pth oc p^ de- 
creasing with increasing (3^; s = 0.96, 0.77, and 0.54 
for Po = 0.2, 2, and 20, respectively, for p e [10-^ lOl. 
The average of /3th versus p can be seen in Figure [iT 
For perfectly spherical contractions, B ex p^^^, since 
B (X R~'^ due to flux conservation, and p ex due 
to mass conservation. For flow perfectly aligned with 
the field, B (x p^, since no amplification of the field can 
take place, and for flow complete ly perpendicular to the 



field, B (X p^ (eg. Kulsrud 2004). Since /3th = c^p/B 



and for an isothermal equation 01 state Cg is constant, one 
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Figure 7. Power spectra for the logarithm of density, P(lnp, k), 
for t = O.ltff, t = 0.3%, and t = 0.6%(sohd, dashed, and dotted 
hnes) for the /5o = 0.2 simulation. Taking the spectrum of the 
logarithm of density allows us to recover the turbulent state even 
at late times, as the collapsing state is filtered out by the loga- 
rithm. Slopes get somewhat fiatter with increasing mean field; the 
/5o = 0.2 simulation has a slope of —1.0, and the = 2 and 20 
simulations having slopes of —1.4 and —1.53, respectively (neither 
shown here) 




10° 10^ 10^ 

^/^min 



Figure 8. Column density power spectra, P(E,/c). Plot style is 
the same as in Figure |2] 

finds 5 = — 1, —1/3, and 1 for perpendicular, isotropic, 
and parallel flow, respectively. This indicates that flow 
is preferentially aligned in all states, even in the weakest 
field simulation, but the alignment increases with mean 
field strength 

The degree to which the velocity and field are aligned 
is shown in Figure [T2j which shows the mass weighted av- 
erage of the magnitude of the cosine of the angle between 
B and v, 

(|cos^|),. = (^)^, (15) 

where the average is done over only material with density 
in bin pi and normalized to the total mass in that bin. 




Figure 9. Thermal-to-magnetic pressure ratio, /3th, vs. p, colored 
by mass fraction, Fm, at t = t = O.ltff (left column) and t = 0.6% 
(right columnd). Top to bottom, /3o = 0.2,0.2 and 20. Increased 
mean field increases the correlation between (3 and p and increasing 
the slope. High density gas has largely /5th ~ 1, showing a tendency 
towards pressure balance in collapsing gas. 

and 512 bins were used. The grey line shows the expecta- 
tion value of I cos6\ for uncorrelated vectors, S/tt ^ 0.64. 
The solid red, green, and blue lines show /3o = 0.2, 2, and 
20, respectively, averaged for several snapshots around 
t = 0.6tff. Light dotted lines show the constituent snap- 
shots, which demonstrates the extremely short timescale 
on which varies at high density. The /3o = 0.2 simula- 
tion has B and v nearly aligned at low density, which is 
consistent with s = 0.96 ^ 1, while the other weaker field 
simulations show less alignment, even a slight tendency 
for B and v to be perpendicular, and correspondingly 
lower values and larger variance in s. 

The left column of Figure [lO] shows the ratio of 
dynamic-to-magnetic pressure for the early snapshot. It 
can be seen from this figure that the typical gas element 
in the low-density turbulent state in the Po = 0.2 run 
is trans- Alfven, as /3dyn = and the gas is evenly 

distributed around /3dyn = 1, with a peak of the PDF at 
/^dyn = 0.54. The other two simulations are more super- 
Alfvenic, with peak /3dyn at 1.2 and 5.6. 

5.2. Energy Ratios in the Collapsing State 

For the high density collapsing gas, /3th and can 
be seen in the right columns of Figures [9] and [lOj The 
transition between turbulent and collapsing states is seen 
quite clearly in Figure [TT] which shows average /3th vs. p. 
In this figure, at Pth ^ 1 the slope in mean /3th flattens to 
almost zero, showing that on average, collapsing gas is in 
thermal-to- magnetic pressure balance. For p G [100, pr], 



10 




logloP 



logloP 



Figure 10. Dynamic pressure ratio ^dyn for t = 0.1% (left col- 
umn) and t = 0.6tff (right column). Top to bottom, /Sq = 0.2,0.2 
and 20. High density gas density gas is dominated by dynamic 
pressure regardless of mean field strength. 




Figure 11. Average thermal-to-magnetic pressure ratio /3th vs. p 
for all three simulations, /3o=0.2, 2, and 20 shown in red, green, and 
blue respectively, at two snapshots, t = O.lfff (solid) and t = 0.6tff 
(dotted.) The change in behavior, from turbulent to collapsing, 
is apparent at /3th ~ 1- The horizontal line shows /3th = 1, and 
vertical lines show the three density regimes, as in Figure [2] This 
highlights the fact that density alone is not enough to descnoe the 
transition from turbulent to collapsing gas. 




Figure 12. Angle between B and v (|cos^|) for a number of 
late-time snapshots (t = 0.57, 0.59, 0.61, 0.63, 0.65fff, light dotted 
lines) and the average (red, green, and blue for /3o = 0.2, /3o = 2, 
/3o = 20, respectively). The horizontal line shows the expectation 
for random 0, while the two vertical lines separate turbulent (left), 
collapsing (center) and unresolved (right) states. 



s = 0.15,0.19, and 0.23 for f3o = 0.2,2, and 20, respec- 
tively. A lower limit of p = 100 was used for these fits 
as this is approximately where the slope in the = 0.2 
simulation flattens dramatically. However, the density 
at which this transition in slope happens is found to de- 
crease with increasing magnetic field, with p ^ 10 and 
0.1 for /3o = 2 and 20, respectively. This indicates that 
density alone is not a sufficient variable to mark the tran- 
sition from turbulent to collapsing gas. 



As shown in the right column of Figure 10, the collaps- 
ing gas is dominated by dynamic pressure for all three 
values of /^q. This is also true for the unresolved gas, 
which while one cannot quantitatively trust these results 
as they are contaminated by numerical resolution prob- 
lems, the likelihood that increased resolution would de- 
crease this ratio by two orders of magnitude is low. This 
suggests that collapsing state is formed from gas that is 
initially super- Alfvenic, where magnetic energy support 
is insufficient to resist ram pressure from the gas, causing 
density peaks that can become gravitationally bound. As 
the mass fraction of gas that is super- Alfvenic decreases 
with this also has the effect of decreasing the amount 
of gas that can collapse to high densities, in turn de- 
creasing the rate and efficiency of star formation. In 
non-self-gravitating turbulence results, the existence of 
super- Alfvenic high-density ^a s is seen even in sub-sonic, 
sub-Alfven simulations (Burkhart et al.|2009), a s well as 



in thermally unstable trans- Alfven simulations (Kri tsukj 
et al.|[2Qllb ). Gravity, however, has the effect of highly 



concentratmg the high density gas in the super- Alfvenic 
regime. This will be compared with observations in Sec- 
tion [S 

6. MAGNETIC FIELD PDF 

Figure [Ts] shows PDF of the magnetic field for t = O.ltff 
(left panel) and both t = O.lt^ and t = 0.6tff (right 
panel). The left plot is linear, and shows V{b)^ where 
b is the fluctuating field, b = \'B — Bo\. The right plot 
is logarithmic, and shows the PDF for the full magnetic 
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field V{B). 

6.1. Magnetic Field PDF in the Turbulent State 

Kritsuk et al. (2012, in preparation) iiave found that 
tiie high field wing of the PDF of magnetic field in 
an isothermal turbulent gas can be well described by a 
stretched exponential, of the form 

V{b)db = cib^'-^/b^Q) exp [-{b/boY] db, (16) 

with a stretching exponent c ~ 1/3. A stretched ex- 
ponential describes a sequence of multiplicative events, 
where 1/c is the depth of the Jiierarchy of events (Frisch & 
|SQrnette||1997t [Laherrere & Sornette 1998). This is not 
unlike the sequence of multiplicative shocks that gen- 
erates the density PDF V{p). For /3o = 0.2 we find 
6o = 3.1 X 10-3 and c = 0.33 for b e [50, 80]; for /3o = 2 
we find bo = 6.1 x 10"^ and c = 0.32, for b G [50, 100]; 
and for /3o = 20, we find bo = 9.0 x 10"^ and c = 0.33 
for b G [50, 120]. Fits to these lines can be seen in black 
along each curve in the left panel of Figure [13] To better 
demonstrate the fit. Figure [l4| shows V{b) against 6^/^, 
restricted to the interval of the fit. Thus we find that 
the number of multiplicative events is the same for all 
three simulations, 1/c ~ 3, but the characteristic scale, 
6o, increases with Po due to the greater ease of inducing 
fiuctuations in a weaker mean field. 

6.2. Magnetic Field PDF in the Collapsing State 

The collapsing state exhibits a power-law V{B) oc B^. 
As an example, for the mid-strength field, /3o = 2 we find 
p = —4 for B G [52, 530]. This fit range was determined 
by the average magnetic field spanned by densities in the 
range p G [pt^ pr]. We compute the exponent qin B (x p^^ 
and find q = 0.43 for this simulation. Combining this 
with m = —1.7 found in Figure [3j we predict p = —3.9, 
in reasonable agreement with the measured value. 

The slope of V{B) increases with increasing mean field 
strength. For the Po = 20 case, we find a fit exponent 
p = -3.22 for B G [27,670], while for po = 0.2 we find 
p = —5.42 for B G [36,451]. Using measured values of 
q and m, we predict slopes for V{B) = —4.5, —4.0, and 
—3.3 for /3o = 0.2, 2, and 20, respectively. The predicted 
value of p for the po = 0.2 run is somewhat lower than 
the measured value. This is likely due to the fact that 
the run had not yet fully developed the collapsing state. 

7. VELOCITY AND ENERGY POWER SPECTRA 

The velocity power spectrum, P('u,/c), is historically 
one of the most studied quantities in turbulence mod- 
eling. Figures [15] and [16] show compensated velocity 
spectra. One curious feature is that the collapsing state 
doesn't show up in the velocity spectra, and only shows 
a mild feature in the kinetic energy spectra, P{p^^'^v^ /c), 
and magnetic energy spectra, P{B^ /c), as seen in Figure 

im 

7.1. Velocity and Energy Power Spectra in the 
Turbulent State 



snapshots (black line). The only variation comes from 
the first O.ltff, where the variation is due entirely to 
the change in the effective spectral bandwidth between 
PPML, the solver used for the initial conditions, and the 
solver in Enzo, which was used for the self- gravitating 
AMR simulations. The rest of the snapshots, nearly in- 
distinguishable from the mean, show almost no evolution 
at all. The other two simulations (not shown) show even 
less variation in the initial phase. 

Figure [16] shows compensated power spectra for all 
three siniulations, with red, green, and blue showing 
Po = 0-2, 2, and 20, respectively. Averages were taken 
over all snapshots, as the variation for /3o = 2 and 
Po = 20 simulations were even smaller than that shown 
in Figure [15] Clearly, the mean field strength has a sig- 
nificant effect on the spectral slope. While the resolution 
of these simulations is too low to make precise measure- 
ments of the slope in the inert ial range, measuring the 
slope is useful to compare scaling among the three simu- 
lations. Slopes for k/k^i^ G [2,30] are —1.46, —1.58, and 
— 1.80 for = 0-2, 2, and 20, respectively. Though there 
is presently no theory to predict the slope as a function 
of mean field in compressible MHD turbulence, the fiat- 
tening of the slope with increasing mean field seen here is 
consistent with what one can infer from the existing the- 
ories for compressible hydro and incompressible MHD. 
The po = 20 slope is consistent with other supersonic 
hydrodynamic simulation of Kritsuk et al. (2007), and is 
the up per end of the slope predicted by Boldyr ev et al^ 
( 2002| ). The slope for the Po = 0.2 simulation is rem- 
iniscent of_th^_^3/2_j^lue_^^ 

model (Iroshnikov 1964 Kraichnan 1965). This reduc- 



tion of slope with ni ean field strength was also seen by 
Kritsuk et al. (2009), who measure —1.94, —1.62, and 
— 1.51 lor Po = 0.2, z, and 20, respectively. This trend of 
increasing slope wit h mean field was continued by |Lemas^ 
ter & Stone (2009), who measured a slope of —1.38 for 



Only the turbulent state is visible in the velocity power 
spectrum, P(v, k). Figure 15 shows the compensated ve- 
locity power spectra, k^^^P(y^k)^ for all snapshots for 
the Po = 20 simulation (grey lines) and the average of all 



Po = 0.U2, which is fiatter and of lower Po than our strong 
field run. The Po = 2 run seems to be in transition be- 
tween the two, with the slope increasing for wavenumbers 
above k/kmin = 10. 

The flattening of the spectra is consistent with the 
magnetic energy coming into equipartition with the ki- 
netic energy, as demonstrated by Figure [Tt] The mag- 
netic energy is about ten times lower thanThe kinetic in 
the Po = 20 simulation for all wavenumbers, so the sim- 
ilarity between this run and other hydrodynamic runs 
is expected. The Po = 2 simulation displays equiparti- 
tion only in a very small band around k/kmin = 10, and 
at higher wavenumbers, where the magnetic energy is 
still below equipartition, the slope of the velocity spec- 
trum seems to increase. The Po = 0.2 simulation has 
near equipartition for the majority of the low to mid 
wavenumbers, and has the flattest spectrum overall. 

Another illuminating feature of Figure [iT] is the sub- 
equipartition nature of the magnetic fleld in these simu- 
lations at small scales at early times (solid curves). This 
lack of equipartition at small scales is somewhat surpris- 
ing, as standard expectation of a small scale dynamo 
in incompressible MHD is to flrst grow exponentially at 
small scales until equipartition is reached, then linearly 
at larger scales, wi th the equipartition wavenumber de - 
creasing with time (Brandenburg fc Subramania n|2005 ). 
However, as the turbulence in this simulation is doni- 
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Figure 13. PDF of magnetic field strength. (Left) PDF of fiuctuating b = |B — So|, V(b), for t = O.lfff. Fits to a stretched exponential 
are shown in black for each line. Note that this figure is linear in b. The smaller volume fractions achieved by the (3o = 20 simulation are 
due to the more rapid collapse, thus earlier refinement, of that simulation. (Right) PDF of the full magnetic field strength, B, with lines 
colored the same as in Figure [2] Here we see a power-law developed at late times due to the collapsing gas. This figure illustrates how 
turbulence and gravity leave dmerent signatures in the PDF of magnetic field strength; at high field strength, a stretched-exponential tail 
is generated by turbulence, while a power-law tail is generated by gravity. 





Figure 14. The stretched exponential section of V(b) at t = O.ltff , 
restricted to the fit range b G [50, 100], and plotted against b^^^. 



Figure 15. Compensated velocity power spectrum k^^^P(v, k) for 
15 snapshots between t = 0.0 — 0.65% for the = 20 simulation. 
All the evolution seen can be attributed to relaxation from the 
initial conditions, no evolution due to gravity is seen. 



inated by shocks, one cannot apply the same physical 
arguments, since the shock jump conditions do not im- 
ply simple equipartition between kinetic and magnetic 
energy in the shock-compressed layer. The increased 
fraction of kinetic energy in compressive motions (see 
Section [8| generates less vorticity per unit of energy 
than sub-sonic turbulence, which in turn gen erates less 
magn e tic ene rgy. This has been explored in Federrath 
et al. (2011a), who demonstrated that as Mach number 
mcreases, the ratio of magnetic to kinetic energy tran- 
sitions from near unity for < 1, to a few percent 
for > 1. This is consistent with our findings here, 
though predominantly at small scales. In our /3o = 2 and 
/3o =0.2 simulations, the imposed large scale field allows 
equipartition to be reached at intermediate scales. 



7.2. Velocity and Energy Power Spectra in the 
Collapsing State 

As seen through the time evolution of P(v, k) in Fig- 
ure [15) the collapse state does not leave a signature on 
the velocity power spectrum. This is due to the volume 
weighted nature of power spectra. The collapsing state 
occupies a very small volume fraction, as seen in Figure 
[2j Because of this, for a signal to appear in the power 
spectrum the values must be extremely large. The val- 
ues of velocity reached by the gas do not vary by the 
many orders of magnitude that the density does in the 
collapsing gas, as seen in Figures [2] and |6] 

If we relate the power spectra fou nd he re with the self- 
similar solutions discussed in Section|3]2) we find that the 
pressure-free collapse gives velocity scaling exponent v = 
+5/7; the expansion wave solution gives = 0, which 
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Figure 16. Velocity power spectra P{v, k) for all three simula- 
tions, averaged over all times. Slopes for A://cmin ^ [2, 30] are —1.46, 
-1.58, and -1.80 for (3o = 0.2,2,20, respectively. 

gives no visible signal to the velocity scaling; finally the 
Larson-Penston solution, with constant velocity, gives no 
contribution to the velocity spectrum. 

While the collapsing state leaves no signature in the 
velocity spectra, it does leave a signature at high k/kmin 
on the kinetic and magnetic energy spectra, P(p^/^'U, /c), 
and P{B^k). In Figure pT] we show three snapshots for 
the energy spectrum, P(p^/^'u, k) (upper set of curves) 
and P{B^k) (lower set), for all three simulations (top 
to bottom, Po = 0.2, Po = 2, /3o = 20). Solid, dashed, 
and dotted lines show t = 0.1, 0.3, 0.6tff , respectively. 
The kinetic spectrum shows an imprint of the collapse in 
the high k/k^i^ > 100 gas, due mostly to the density- 
weighted nature of this statistic, and the fact that, on 
small scales, the density increases by as much as 5 orders 
of magnitude. Strictly speaking, the magnetic spectrum 
is also a volume weighted quantity, so the increase in 
magnetic energy seen is a product of additional field am- 
plification from the collapse itself. However, as discussed 
in Section [5| B oc p^, with q ^ 0.4, which is a similar 
to the density dependence of kinetic energy, so P{B^ k) 
can be considered to be implicitly density weighted. The 
initial (turbulent) field distribution is a product of the 
shock dominated turbulence, as discussed in the previ- 
ous section. After gravity is introduced, the field sees 
additional amplification due to the collapse itself. This 
field amplification is again consistent with increased com- 
pressibility due to increased /3o, as the increase in P{B) 
at high k/kmin is greater for the weaker field runs. 

If one examines the energy scaling in the context of 
self-similar spheres, one finds that P{p^^'^u) scales as k~^^ 
for the Larson-Penston and pressure-free solutions, and 
for the expansion wave solutions. These increasing 
solutions are not seen in Figure [T7[ Again this is likely 
due to the fact that the superposition of the turbulent 
and collapsing spectra makes distinguishing between the 
two difficult. 

8. HELMHOLTZ DECOMPOSISITION 

The formation of high density material necessarily 
comes from compressive motions. A standard method 
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Figure 17. Energy power spectra, P{p^^'^v^k) (upper curves in 
each panel, red, green, and blue) and P{B,k) (lower curves in each 
panel, cyan, magenta and yellow) for /3o = 0-2 (top panel) ^^o = 2 
(middle panel) and /3o = 20 (bottom panel). Three different times 
are shown; t = O.ltff (solid line), t = 0.3tff (dashed line), and 
t = 0.6% (dotted line) 
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of examining the compressional and solenoidal modes of 
the velocity field is the Helmholtz decomposition. We 
split the velocity field, v, into two components 

(17) 
(18) 
(19) 



V = Vc + Vs, 

V • V, = 0, 

V X Vc = 0, 
where we find Vc and Vg in Fourier space, 

Vc = (fc • \-)k, 

Vs = V - Vc 



(20) 
(21) 



where v is the Fourier component of v, and k is the unit 
wave vector. We then examine the ratio of the power 
spectra x(^) = P{^c-, k)/P{'Vs^ /c), which gives us a scale- 
wise measure of the importance of compressibility of the 
gas. 

8.1. Helmholtz Decomposition in the Turbulent State 

The ratio of compressible to solenoidal power, xl^)? 
is shown in Figure [l8] for the same fiducial snapshots 
as in Figure [2] Forpo = 20, we find x(^) ^ 0-4 in 
the range k/W^in G [2,30], which is less than 0.5 as 
expected from purely geometrical considerations: for a 
given wave vector k there is only one longitudinal (com- 
pressional) component, while there are two transverse 
(rotationa l) compo nents. This is consistent with the find- 
ings of Kritsuk et al.| (12010b), where the increase of the 
magnetic field strength decreases x(^)- Again we find 
the increased compressibility with increased Pq. At low 
k/kmin < 50 there is a monotonic increase in x(^) with 
/3o. The f3o = 20 simulation also shows a substantial low 
k/kmin < 2 compression at late times, seen to a lesser 
extent in the other two. The Po = 2 run shows a peak 
in x(^) at k/k^i^ ~ 10, probably with a similar ori- 
gin as the k/kmin < 2 increases, namely a result of the 
large scale forcing pattern. In all cases, there is a gen- 
eral trend of increasing compressibility at higher k/kmin- 
As discussed in Section [7[ this increase in compressibil- 
ity causes reduced amplmcation of the magnetic field at 
high k/kmin- 

8.2. Helmholtz Decomposition in the Collapsing State 

There is no substantial increase in compressive mo- 
tions at high k/kmin as time progresses, indicating, as for 
the velocity power spectrum, P{v,k), that this statistic 
is not especially sensitive to gravity, due to its volume- 
weighted nature. It seems that the lower levels of mag- 
netization in the Po = 20 simulation is enough for that 
run to have some increased compressional motion at high 
k/kmin- The stronger field runs show almost no evolu- 
tion. At low k/kmin^ all three simulations show an in- 
crease in power at k/kmin < 2. This is possibly due to 
the effects of gravity causing large scale contraction of 
the gas. 

9. DISCUSSION 

The combined effects of gravity, turbulence, and mag- 
netic fields have yet to be incorporated in star formation 
theory with their proper respective weights. Here we dis- 
cuss mo dific ations to two recent star formation models 
(Section [9?T] ) and interpret several recent observations in 
this light (Section [9^. 




^/^min 



Figure 18. The ratio of compressible to solenoidal motions, 
x(k) = P(vc, /c)/P(vs, /c), for all three simulations (red, green 
and blue for /3o = 0.2, /3o = 2, = 20) and two snapshots 
{t = 0.1, 0.6%, solid and dotted lines, respectively) 



9.1. Implications on Theory 

In the past decade, the lognormal density PDF of tur- 
bulence has been used to predict a number of properties 
of star formation, including the star formation rate. Here 
we will examine the implication on two works in part icu- 
lar, n amely the star formation rate models of Kru mholz 
'ITMcKee (2005), hereafter KM05, and Pad oan fc Nord- 
llund (2011), hereafter PNll. In both of these works, the 
star formation rate per free fall time, SFRff , is predicted 
to be proportional to the mass fraction above some crit- 
ical density, using a lognormal distribution. Thus 



SFRff = e / pV{p)dp 

'per 

"(T^ — 21npc 



e 1 + erf 



(22) 
(23) 



where we have used Equation 11 in the derivation of the 
second of these. This model presumes that the effects 
of gravity do not alter the structure of the star-forming 
cloud. However if the effects of gravity are present dur- 
ing the formation of the cloud, it is not unreasonable to 
suppose that a power-law density PDF might be present 
before the onset of turbulence. Here we will explore what 
effects the distribution may have on the star formation 
rate. 

KM05 and PNll differ in two ways: the nature of e 
and the nature of per- The parameter e contains the 
timescale for collapse, as well as the fraction of collaps- 
ing gas reaching the proto-star. The other difference is 
the nature of pcr- In KM05, per is the density at which ve- 
locity perturbations become subsonic and lose their tur- 
bulent support. In PNll, pcr is the post shock density 
that yields clumps that are larger than the sum of the 
Bonnor-Ebert and magnetic critical masses. 

As we have shown in Section [7| the scaling of the tur- 
bulence depends on the mean magnetic field. This in 
turn has consequences for the expected critical density. 
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The KM05 result predicts that 

A, = Lo(Cs/<To)l/^ 

thus 



(24) 
(25) 



M'^^ (26) 

^0 / 



where Ag is the sonic length, at which velocity pertur- 
bations become subsonic, Lq is the cloud size, ctq is the 
velocity fluctuation at the outer scale, and (/) is a nu- 
merical factor of order unity. In their work they find 
(j) = 1.12. We can relate q to the velocity spectral index 
uhy q = — (z/+l)/2. For z/ = — 2, as expected for pressure 
free turbulence, q = 1/2, while for u = —1.5 as seen in 
our Po — 0.2 run, q = 1/4. If we take A4 = 10, we find an 
increase of several orders of magnitude in per when using 
the spectral scaling of the /3o = 0.2 simulation over the 
f3o = 20 simulation, assuming everything else stays con- 
stant. In PNll, Per does not depend directly on z/, but 
will likely impact the post-shock magnetic field distri- 
bution, similarly increasing p^r for increasing mean field 
strength. 

Another impact of our results here is the result of re- 
placing the turbulent PDF, V{p), with one that transi- 
tions to a power-law at some transition density pt. Thus, 



Vp{p)dlnp ■■ 



N- 



[Npop'^dlnp^, 



Qnp- 



dlnp, p<pt 

P > Pt, 
(27) 

rather than a plain lognormal. Here the normaliza- 
tion N is determined by the normalization requirement. 



= 2 1 + erf 



1, and is given by 



23/2^ 



(28) 



From here we are left to determine the critical density for 
star formation, pcr, and the density at which the PDF 
transitions from lognormal to power-law, pt. These are 
not necessarily equal; the forces of rotation and mag- 
netic fields will suppress the collapse of material once 
the gas is within the potential well of the self-similar 
sphere, whereby Pcr may be larger than pf If we assume 
that V{p) is continuous and different iable, we can piece 
together an analytic estimate for po and pt. Differentia- 
bility is not necessarily the correct condition to take, but 
in the absence of a better condition for the transition 
density pt, an assumption must be made. The PDF in 
Figure [2] as well as some observations do appear to be 
different iable. Nevertheless, we can use this assumption 
to examine the sensitivity of SFRff to pt and po. These 
two conditions give us 



pt = -{2m-l)a^ 



(29) 

^g.^m(m-l)/2^ (30) 

Using the fit values from Table [T] we find values listed 
in Table [3) We also compute values for the Larson- 
Penston sphere, with p oc and a = 1.5 as found 



for the majority of our simulations, and a pressure- free 
collapse with a taken from Equation 12 for b = 0.4 and 
M = 9. 

The increase in predicted star formation rate is found 
by using this piecewise PDF instead of a lognormal in 
Equation 23 In Figure [19) we plot the piecewise cu- 



mulative mass fraction in Equation [27] relative to the 
purely lognormal mass fraction versus critical density Pcr- 
We have used parameters from all three simulations, red, 
green, and blue lines for /3o = 0.2, 2, and 20, respectively, 
the Larson-Penston solution in black, and the pressure 
free solution in grey. We find that the increase in star 
formation rate can be quite large, and is quite sensitive 
to the details of pt and m. This clearly predicts the in- 
correct behavior in rate relative to /3o, as the /3o = 0.2 
model has a measured collapse rate lower than the other 
two (see Figure [5| . This is due to the selection of pt, 
which is significantly lower for Po = 0.2 than the other 
two, caused by the lower value of a. While this is an 
incorrect prediction, it demonstrates how sensitive this 
type of model is to the parameters of the fit. In reality, 
per and pt will be functions of a stability criterion that 
fully incorporates magnetic fields and turbulence . 

This result is consistent with the findings of Ch o fc Kim| 
(2011 ), who used 512^ unigrid simulations to examine the 
collapse rate for various values of critical density. For the 
largest critical values, pcr = 500po5 they find an increase 
of 2400 relative to the predicted value of KM05. 

Future work will analyze bound clumps to determine 
a proper value for p^r and pt, and further evaluate this 
prediction as well as KM05 and PNll. 



Table 3 

Parameters for extended density 
PDF 



model 


m 


<j 


Po 


Pt 


0.2 


1.2 


1.75 


2.4 


6 


2 


1.5 


1.75 


3.2 


17 


20 


1.5 


1.64 


2.4 


13 


PF 


1.8 


1.71 


3.4 


24 


LP 


1.5 


1.5 


1.9 


10 



One question that is not addressed by this work is when 
star formation begins in the lifetime of a cloud. That is, 
are the effects of self gravity already present before the 
lognormal PDF is established, in which case a power-law 
V{p) is the best picture to take, or is a star- forming cloud 
turbulent first, and then transitions to_self-gravit at ing 
through cooling? Kainulainen et al.l (|2009]) two families 
of clouds, one with only lognormal PDFs and no star 
formation, and one family that includes power-law tails 
and high rates of star formation. However it is presently 
unknown if the star formation begins before or after the 
power-law tail, or if the tail impacts the star formation 
itself. 

9.2. Interpretation of Observations 

The role of magnetic fields has been vigorously debated 
over the last few decades. Initially magnetic fields were 
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tion in field strength with scale. Crutcher et al.l (I2 Q1Q[ ) 

ti tor 



Figure 19. Excess cumulative mass fraction of a power-law PDF 
relative to a lognormal PDF. This gives the amount by which the 
two primary star formation rate models are low. The black line 
uses a = 1.5, m = — 1.5,po = 0-5. Red, green, and blue lines use 
(J taken from Table ^ m from Figure [3] and pQ from Equation |30| 
The grey line uses a pressure-free solution (PF) 



the dominant mechanism for regula ting star formation 
( Mouschovias 1976 Shu et al. 1987), then relegated to 
an m signitica nt role as supersonic turbulence took hold 
^lac Low & Klessen 2QQ4[ [Krumholz fc M cKee 2005). 
Recent observations have given a mixed review of the role 
of magnetic fields, with some measurements indicating 
a strong field, and some indicating a weak field. Here 
we will discuss the implications of our present model on 
re cent obse rvations . 
Li et al.| (|2009) measured starlight polarization in a 



number of molecular clouds, and averaged said polar- 
ization over large (30 pc) scales of the cloud and small 
(0.3 pc) scales centered on high density cloud cores. 
The resultant alignment of magnetic fields on large and 
small scales is interpreted to imply sub-Alfvenic magnetic 
fields, by way of comparison to a set of super-Alfvenic 
(A^A = 2) and sub-Alfvenic {M.p^ = 0.7) simulations. In 
sub-Alfvenicturbulence, the kinetic energy of the gas is 
unable to alter the alignment of the magnetic field at all 
scales, leading to alignment with the mean field at all 
scales and a correlation between direction at both large 
and small scales. In contrast, super-Alfvenic turbulence 
can greatly alter the magnetic field direction, reducing 
the correlation between large and small scales. This is 
consistent with the increase in mean alignment between 
V and B with decreasing (Figure fl2|, as in our trans- 
Alfvenic simulation there is a trenoTowards alignment 
that is absent from the higher /3o runs. In our trans- 
Alfvenic simulation, kinetic energy is insufficnent to alter 
the direction of the field relative to the mean, and in turn 
velocity and magnetic vectors are aligned. In the super- 
Alfvenic simulations, on the other hand, kinetic motions 
dominate and the field direction can be altered by the 
velocity. On the other hand, recent measurements of the 
directions of outfiow from T-Tauri stars are uncorrelated 
with mean magn etic field ( Menard fc Duchene|2004 Tar- 



gon et al. 2011), though there is a weak correlation fo r 
earlier Class and Class 1 objects (Targon et al. [2011 ). 
Finally, two sets of observational papers show a transi- 



find that density is uncorrelated with field strength for 
densities nH < 10^, and correlated as B oc n^-^^ for 
higher densities, using Zeeman splitting. |Heyer fc Brunt 



(12012) use velocity anisotropy in the Taurus molecular 
cloud to show that the low column density gas, with more 
anisotropic fiows, is likely trans- Alfvenic, while high col- 
umn density gas, with more isotropic fiow, is likely super- 
Alfvenic. 

This suite of measurements is consistent with what we 
have seen in this study, and can be understood by the 
combination of Figures [lO| [T2| and [T7| focusing on the 
/3o = 0.2 simulation. Through the action of supersonic 
turbulence, the effect of magnetic fields is a function of 
scale. Low k/kmin gas, where the turbulence is gener- 
ated in the presence of an appreciable mean field, comes 
into equipartition. This allows low density gas to ex- 
hibit alignment seen in observations. However in mate- 
rial at higher k/kmin^ where the structures are generated 
by shocks, the field is unable to come into equipartition. 
This high density gas is further selected to be super- 
Alfvenic in post-shock gas where the magnetic energy is 
too weak to resist compression. 

We also find that our simulations are filled with fila- 



mentary structures, as seen by Men'shchikov et al. (2010) 
with Hershel in the Aquila and Polaris clouds, see b'ig- 
ures [l] and [20] Figure [20] shows a restricted region of 
the t = 0.6tff snapshot ofthe /3o = 2 simulation, over- 
laid with magnetic field vectors (right panel) and velocity 
vectors (left panel). High density material, with p > 10, 
is nearl y all concentrated in filam entary structures, as 



Arzoumanian et al. 



2011). Further analysis is 



necessary to quantify the filament width distribution 
as in that work filaments cluster tightly around O.lpc. 
Magnetic fields are weak and perpendicular to the struc- 
ture with the cores, while in the lower density filamen- 
tary tail the field is parallel to the structure. This is 
expected by the fact that all high density material is 
super-Alfvenic: the filament aligned with the magnetic 
field will have maximally amplified the magnetic field, 
halting the collapse, while the material with B and v 
more closely aligned has lower supporting pressure due 
to the reduced field amplification, and can go on to form 
"stars" . 

We can make further contact between the collapsing 
state and observations of prestellar cores seen in the 
Aquila star-forming cloud, as reported by And re et al. 
(2010). We can estimate the mass-size relation from a 
self- similar sphere as 



M{R) 



= 47r / 
Jo 



p{r)r^dr ex K 



in+3 



(31) 



which gives us M(r) ex for the LP solution, M(r) (x 
'^^•^ for the PF solution (see also [Kritsuk fc Norman 



(2011)). The value found in Aquila is 1.13 (Philippe 
Andre, private communication.), which is somewhat 
shallower than that of the LP solution. For complete- 
ness, the expansion wave solution predicts M(r) oc r^-^, 
though is less relevant due to the fact that the expansion 
wave solution presupposes a central singularity, while the 
Andre result focused on starless cores. 

10. CONCLUSIONS 
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Figure 20. Projections of a sub set of dense gas from the (3o = 2 simulation at t = 0.6tff . Overlaid are magnetic vectors (left) and velocity 
vectors (right). The dense parts of this filament lie at the convergence of opposing velocity streams, where the magnetic field is roughly 
aligned to the velocity. Higher in the stream, where the field is perpendicular to the velocity, collapse is suppressed. 



In this work, we use high resolution AMR MHD sim- 
ulations to examine the combined effects of self-gravity 
and magnetic fields on supersonic turbulence. We find 
that self-gravity bifurcates the cloud into two phases: 
a low density, turbulent cloud; and high density, self- 
gravitating cores. The two phases have substantially dif- 
ferent statistical properties. The magnetic field serves 
to effectively decrease compressibility of the gas as the 
mean magnetic field is increased. 

The simulations were perfor med using th e AMR Enzo 
code (Bryan et al. 1995; O'She a et al.|2QQ4[ ), extended to 
MHD by |Collins et al., ( ,2QlQ^ . 'I'he three simulations use 
an rms Mach number JVl = 9, virial parameter a = 1, 
and three values of the mean magnetic field, giving initial 
plasma beta = 0-2, 2, 20. The coarse level had 512^ 
zones, and 4 levels of refinement were used, refining to 
keep the Jeans length resolved by 16 zones. This gives 
an effective resolution of 8196^. 

The low density turbulent state exhibits properties an- 
ticipated by earlier works in supersonic hydro and MHD 
turbulence The effects of increasing mean magnetic field 
strength are to decrease both the compressibility and the 
ability of the gas fiow to bend field lines. The properties 
of the turbulent state can be summarized as follows: 

• Density PDF, V{p)^ can be roughly represented by 
a lognormal, given by Equation [ll] The width of 
the lognormal increases somewhat with increasing 
mean thermal-to-magnetic pressure ratio, /^q. 

• The density power spectrum, P(p, /c) ex , tends 
to fiatten with increasing Mach number. We find 
slopes that are fiatter than other measurements of ( 
in supersonic turbulent simulations, in proportion 
to the higher Mach number relative to earlier work. 
We find that ( weakly decreases with increasing 
Pq. For the early snapshot, ( = —0.42, —0.58, and 
—0.62 for /3o = 0.2, 2, and 20, respectively; 

• The column density power spectrum, P(i;, k) ex 
k^, is also fiatter than observed. We find A = 



-1.39,-1.61, and -1.53 for (3o = 0.2, 2, and 20, 
respectively, while values of A < —2.25 have been 
reported elsewhere. This discrepancy is likely due 
to the decreased density range in earlier measure- 
ments, in part due to the reduced dynamic range 
available to the single tracer molecules used for the 
observations. 

• Thermal-to-magnetic pressure ratio, /3th, shows a 
decreasing correlation between /3th and p with in- 
creasing f3o. The average slope also decreases with 
decreasing field. For /3th cx p^, s = 0.97, 0.77, and 
0.54 for f3o = 0.2, 2, and 20, respectively. 

• Dynamic-to- magnetic pressure ratio, /3dyn, shows 
only mild decrease in correlation between /3dyn and 
p, and the peak of the distribution decreases with 
increasing Pq. 

• The magnetic field and velocity are nearly aligned 
for mid- to low-density gas in the strong field case, 
and decorrelated for the other two cases. 

• the PDF of magnetic field, V{b) cx exp[-(6/6o)1, 
shows a stretched exponential with a stretching ex- 
ponent of c ^ 1/3, as given in Equation [T6j for the 
high field strength gas. The slope decreases with 
/3o, with bo = 0.15,0.19 and 0.23 for /3o = 0.2, 2, 
and 20, respectively, and the peak decreases with 
increasing /3o. 

• The velocity power spectra, P{v^k) cx /c^, show 
an increasing slope as /3o increases, with the weak 
field case resembling supersonic hydrodynamic tur- 
bulence; we find that u = —1.46, —1.58, and —1.80 
for /3o = 0.2, 2, and 20, respectively. 

• Equipartition between kinetic and magnetic energy 
is only seen at large scale for the strongest field 
case, /3o = 0.2, and a small range of intermediate 
wavenumbers for /3o = 2; the magnetic energy is an 
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order of magnitude or more for all high wavenum- 
bers in all simulations, and all wavenumbers in the 
weak-field case, = 20. 

• The ratio of compressive-to-solenoidal motions, 
x(A:), increases with Pq^ indicating increased com- 
pressibility with decreased magnetic field strength. 

The high-density collapse phase exhibits properties 
consistent with spherically symmetric isothermal collapse 
solutions, with power-law density profiles p ocr^. Again 
the effect of the magnetic field is to decrease the com- 
pressibility of the gas. The properties of the collapsing 
state can be summarized as follows: 

• The density PDF is well approximated by a power- 
law, V{p) (X p^, where m = 3/n, and m = 
-1.80, -1.78, and -1.65 for Po = 0.2, 2, and 20, 
respectively. The mass flux rate, which can be seen 
as a proxy for star formation rate, is an increasing 
function of /Sq. The slope in the strong field case is 
continuing to increase at the end of the simulation, 
due to the slower collapse in this case. 

• The density power spectrum, k) oc k^^ shows a 
positive slope, consistent with the expected behav- 
ior of a self-similar solution, ( = — 2(m + l), though 
detailed comparison to theory is difficult due to the 
superposition of the turbulent state. We find that 
C = 0.86,1.12, and 1.2 for /3o = 0.2, 2, and 20, 
respectively. 

• The column density spectrum P(S, k) (x k^^ shows 
flat spectra, with A ~ for all simulations; this 
behavior is potentially observable, provided obser- 
vational tracers with high enough dynamic range 
are used. 

• The thermal-to-magnetic pressure ratio, ^th, is 
near unity for all values of Pq. 

• The dynamic-to-magnetic pressure ratio, Pdyn^ is 
larger than unity by at least two orders of magni- 
tude for all values of ^^dyn- 

• The high density gas shows a decrease in the mean 
alignment between B and v with increasing den- 
sity. 

• The magnetic field PDF shows a power-law behav- 
ior, V{B) oc BP, with p = -5.42, -4, and -3.22 
for /3o = 0.2, 2, and 20, respectively. 

• The velocity power spectra show no evidence of the 
collapsing state, due to the small volume fraction 
of the collapsing gas and lack of severe increase in 
velocity, as seen in the density power spectrum. 

• The kinetic energy spectra, P(p^/^iA, /c), and 
P{B,k) show increases in power at high /c/Zcmin, 
consistent with the density weighted nature of these 
two statistics. 

• The ratio of compressive-to-solenoidal power sim- 
ilarly shows almost no signature of the collapsing 
state. 



We find that for certain values of the transition from 
lognormal to power-law density PDF, the feedback on the 
predicted star formation rate can be quite large. This 
has sensitive dependence on the transition density, pt 
between lognormal to power-law, and the critical den- 
sity, Per, above which cores are considered to be gravi- 
tationally bound. Additionally the fiattening of the ve- 
locity power spectrum in the trans- Alfvenic simulation 
will cause high density material to remain gravitation- 
ally bound to smaller scales, increasing pcr from values 
expected from hydrodynamic scaling. 

High density collapsing gas is created very super- 
Alfvenic, with A4a ^ 100, even in the trans- Alfvenic sim- 
ulation. In the trans-Alfvenic simulation, we find large 
scale imprints of the magnetic field, and at small scale 
this imprint weakens. This is consistent with a number 
of recent observations of the properties of star-forming 
clouds. High density cores seem to be found primarily in 
filaments, in regions where the magnetic field is aligned 
perpendicular to the filament. Filaments with longitu- 
dinal magnetic fields do not show high density material, 
consistent with the super- Alfvenic nature of the collaps- 
ing state being a necessary condition for collapse. 
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